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The  development  of  an  air-breathing  hypersonic  vehicle  employing  scramjet  propulsion 
is  an  ongoing  research  endeavor.  Because  of  high  velocity  (>Mach  5),  length  and  positioning 
of  the  engine,  and  relative  sleekness,  the  flexibility  of  the  vehicle  is  significant  and  there 
are  strong  couplings  between  thrust  and  pitch  .  As  a  result,  control  design  for  such  a 
vehicle  is  a  challenge.  In  previous  works,  linear  controllers  have  been  designed  for  a  model 
of  the  longitudinal  dynamics  of  a  specific  air-breathing  vehicle  possessing  the  same  number 
of  inputs  and  outputs.  In  this  paper  we  consider  a  control  design  for  the  same  vehicle 
model,  but  we  restrict  our  attention  to  controlling  only  two  outputs,  namely  the  altitude 
and  velocity,  while  we  employ  as  control  inputs,  the  elevator  deflection,  total  temperature 
change  across  the  combustor  and  the  diffuser  area  ratio  of  the  combustor.  The  specific 
control  problem  addressed  in  the  paper  is  the  design  of  a  controller  that  ensures  asymptotic 
tracking  of  altitude  and  velocity  reference  trajectories,  while  using  the  redundancy  in  the 
inputs  to  optimize  the  performance  in  steady-state.  As  a  matter  of  fact,  since  the  system 
is  not  square,  the  steady  state  solutions  that  enforce  perfect  tracking  are  nonunique.  The 
controller  employs  a  parameterization  of  all  possible  steady  state  trajectories  that  is  used 
for  optimization  of  the  steady  state  input  while  providing  perfect  tracking  and  fulfilment 
of  constraint  on  the  magnitude  of  the  control  input.  Simulations  results  are  provided  to 
validate  the  proposed  approach. 


I.  Introduction 

Hypersonic  air-breathing  vehicles  have  been  proposed  as  a  promising  technology  for  reducing  the  cost  of 
launching  small  satellites  or  other  vehicles  into  low  earth  orbit.  Furthermore  the  technology  may  present  the 
Air  Force  with  efficient  aircrafts  capable  of  rapid  global  response.  Successful  progress  has  been  made  with 
hypersonic  aircrafts  such  as  NASA’s  X-43A  scramjet  research  aircraft,  which  has  set  new  Guinness  world 
records  for  jet-powered  aircraft  speed  three  times  in  a  row,  reaching  speeds  up  to  Mach  9.6. 1 

Modeling  of  hypersonic  air-breathing  vehicles  is  a  subject  of  active  ongoing  research.2  4  Most  efforts  have 
been  concentrated  on  modeling  only  the  longitudinal  dynamics  since  they  present  a  number  of  challenges  by 
themselves,  especially  due  to  strong  coupling  between  propulsive  and  aerodynamic  effects.  Since  modeling 
is  still  in  its  early  stages,  there  are  very  few  papers  that  have  been  published  on  the  guidance  and  control  of 
such  aircrafts.  Aside  from  the  work  by  Schmidt  and  Velapoldi5  and  Groves  et  al.6  which  focuses  specifically 
on  a  air-breathing  vehicle  model,  the  largest  body  of  work  done  in  this  area  regards  the  control  of  a  generic 
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Figure  1.  Vehicle  Geometry  for  Hypersonic  Air-breathing  Vehicle  from  Bolender  and  Doman2 


hypersonic  vehicle  models  (for  instance,  see  the  robust  nonlinear  controllers  proposed  by  Marrison,  et  al.7 
and  the  adaptive  sliding  mode  controller  by  Xu,  et  al.8). 

As  can  be  expected,  controlling  air-breathing  hypersonic  vehicles  is  a  difficult  task.  Even  after  lineariza¬ 
tion  of  the  dynamics,  control  challenges  are  still  present  such  as  instability,  weak  controllability  of  some 
modes,  considerable  coupling  between  modes,  and  potential  over  actuation  depending  on  input  and  output 
selection. 

We  will  focus  on  control  system  design  for  the  model  developed  by  Bolender  and  Doman2  which  has  been 
previously  used  in  Ref. 6.  The  plant  model  has  three  inputs,  namely  the  total  temperature  change  across  the 
combustor,  elevator  deflection,  and  diffuser  area  ratio.  As  opposed  to  the  case  considered  in  Ref.6,  where  the 
outputs  to  the  be  controlled  were  the  vehicle  speed,  angle  of  attack  and  altitude,  we  consider  in  this  paper 
the  problem  of  controlling  speed  and  altitude  only.  Thus  in  the  case  considered  here  we  have  an  overactuated 
system  and  input  trajectories  that  provide  perfect  tracking  are  non-unique. 

When  considering  an  overactuated  system  it  is  implied  that  for  output  tracking  there  is  an  input  redun¬ 
dancy.  This  redundancy  can  be  used  for  optimization  and  potentially  for  satisfying  constraints.  Traditionally 
this  redundancy  has  been  dealt  with  using  optimal  control  or  control  allocation.9  13  In  this  paper,  we  propose 
a  different  approach,  which  has  been  inspired  by  the  work  of  Johansen  and  Sbarbaro.14  Using  an  exosystem 
to  generate  the  reference  trajectories  to  be  tracked,  the  tracking  control  problem  is  cast  as  a  regulation 
problem.  By  parameterizing  all  solutions  to  the  regulator  equations,15  we  separate  the  control  design  into 
a  design  of  a  steady  state  controller  and  a  stabilizing  controller.  Those  two  controllers  can  be  considered 
separately,  e.g.  using  state  feedback  and  optimizing  the  steady  state  behavior  of  the  system  parameterizing 
the  steady  state  trajectory.  Furthermore,  we  can  exploit  this  parameterization  to  keep  the  input  within 
given  constraints  at  steady  state. 

The  paper  is  organized  as  follows:  Section  II  introduces  the  model  of  the  air-breathing  hypersonic  vehicle. 
In  Section  III  the  tracking  problem  is  introduced  and  all  solutions  of  the  Francis  equation  are  parameterized 
for  an  overactuated  system,  thus  providing  separation  of  steady  state  and  stabilizing  controllers.  Sections 
IV,  presents  an  optimizing  steady  state  controller,  and  in  Section  V  constraints  are  applied  to  the  steady 
state.  Finally,  in  Section  VI  simulation  results  are  presented,  and  conclusions  are  discussed  in  Section  VII. 

II.  The  Model  of  the  air-breathing  Hypersonic  Vehicle 

The  model  of  longitudinal  dynamics  of  an  air-breathing  hypersonic  vehicle  we  consider  in  this  paper  has 
been  developed  by  Bolender  and  Doman2  and  used  for  control  design  in  Ref.6  .  Compressible  flow  theory 
was  used  to  determine  pressure  and  shock  angles,  a  model  of  a  two  cantilever  beam  was  used  to  model  the 
flexible  modes,  and  Lagrange’s  equations  equations  were  used  to  derive  the  equations  of  motion.  Figure  1 
shows  the  vehicle  geometry  used  for  the  model.  The  following  set  of  equations  are  the  non-linear  equations 
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of  motion, 
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where  Vt,  a,  Q,  h,  0  and  rji  are  the  speed,  angle-of-attack,  pitch  rate,  altitude,  pitch  angle  and  ith  generalized 
elastic  coordinates  respectively,  Ft ,  Fl,  Td,  Ah,  and  Ah  are  the  trust,  lift,  drag,  and  ith  generalized  elastic 
forces  respectively,  M  is  the  pitching  moment  about  the  y- axis,  Iyy  is  the  moment  of  inertia,  and  Q,  uji  are 
the  ith  damping  coefficient  and  natural  frequency  of  the  elastic  modes.  Also, 

kl  =  1  +  £’  V’1  =  /°Ll  ^ 

k2  =  1  +  P~,  ^2  =  f°_Lm2€(t>2(e)de 

lyy  *->Z 

where  mi,  m2  are  the  mass  densities  of  the  forebody  and  aftbody  respectively,  and  </>  1,  0 2  are  the  mode 

shapes  for  the  forebody  and  aftbody  respectively.  The  states  and  inputs  are  given  in  Table  1.  The  control 

inputs  are  the  total  temperature  change  across  the  combustor  ATo,  elevator  deflection  5,  and  diffuser  area 
ratio  The  inputs  affect  the  forces  and  moments  which  in  turn  affect  the  states.  As  expected,  AXq  and 
Ad  greatly  affect  the  velocity  and  S  mainly  affects  the  angle-of-attack  and  pitch.  However,  there  exist  strong 
couplings  that  can  not  be  ignored  which  are  caused  largely  by  the  high  operating  velocity  of  the  vehicle,  the 
positioning  of  the  engine,  and  the  flexible  effects. 

Among  the  possible  output  choices,  in  this  paper  we  use  the  vehicle  altitude  and  speed  as  outputs  to  be 
regulated,  i.e.  we  wish  to  have  the  speed  and  altitude  track  a  desired  reference  trajectory. 

The  full  non-linear  model  is  very  complex  and  thus  does  not  lend  itself  to  a  model  based  control  design. 
For  our  design  we  will  therefore  consider  a  linearized  model  which  was  found  using  numerical  methods  to 
linearize  the  non-linear  model  around  a  trim  condition,  which  is  the  one  given  in  Table  2.  The  linearized 
model  is  expressed  in  the  usual  form 


Xp  —  ApXp  H-  BpUp  (9) 

Up  =  CpXp ,  (10) 

with  state  xp  G  M9,  input  up  G  M3  and  output  yp  Gl2.  For  that  specific  trim  condition,  the  linearized  model 


States 

Inputs 

vs 

Vehicle  speed  -  Output 

<5e 

Control  surface  deflection 

a 

Angle-of-attack 

ATo 

Total  temperature  change  across  combustor 

Q 

Pitch  rate 

Ad 

Diffuser  area  ratio 

h 

Altitude  -  Output 

Cowl  lip  position  (Fixed  at  trim) 

9 

Pitch  angle 

Vi 

Generalized  elastic  coordinates 

Vi 

Generalized  elastic  coordinates  time  derivatives 

Table  1.  Definitions  of  the  states,  outputs  and  inputs  for  the  hypersonic  air-breathing  vehicle  modeled  in  Ref.2 
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States  at  Trim: 

Inputs  at  Trim: 

Vi  =  7846.4  ft/sec  771  =  1.6105 

Se=  7.8589° 

a  =  1.9605°  ?)i  =  0 

ATo=  484.49  degR 

Q  =  0  deg/sec  772  =  1.4582 

Ad=  0.35 

h  =  85000  ft  f)2  =  0 

2^=7.0967  fixed 

6  =  1.9605° 

Table  2.  Trim  condition  for  the  inputs  and  states  of  the  nonlinear  hypersonic  air-breathing  vehicle  model. 


has  no  transmission  zeros  and  the  eigenvalues  of  the  A  matrix,  i.e.  the  poles  of  the  plant,  are 

[-0.392  zb  20.1i,  -0.320  ±  16.4i,  0.811,  -0.881,  -1.91  *  10-4,  -4.63  *  10-5  ±  0.0399i]. 


In  order  to  fulfill  actuator  bandwidth  constraints  we  append  to  the  plant  a  simple  model  of  the  actuator 
dynamics,  given  by 


xs  =  Asxs  +  Bsus , 


where 
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The  resulting  augmented  system  is  written  as 


(ii) 


x  =  Ax  +  Bu 
V  =  Cx, 


(12) 


where 


The  linear  system  with  states  x,  input  u  and  output  y  will  be  referred  to  as  the  linear  model  of  the  air- 
breathing  hypersonic  vehicle.  This  is  the  model  that  we  will  use  for  our  control  design  in  the  sequal. 


III.  The  Tracking  Problem 


The  linear  regulator  problem  is  a  fundamental  problem  in  modern  control  which  models  a  variety  of  control 
applications.  Specifically,  the  output  regulation  problem  consists  in  letting  the  output  of  a  system  track  a 
reference  trajectory  or  rejecting  a  disturbance  generated  by  an  autonomous  LTI  system.  The  conditions  for 
existence  of  a  solution  of  the  regulator  problem  can  be  summarized  by  two  conditions:  stabilizability  of  the 
system  and  solvability  of  the  so-called  Francis  equations.15-17 

To  set  up  the  tracking  problem  that  will  be  investigated  in  this  paper,  let  us  consider  the  linear  system 

x  =  Ax  +  Bu 
V  =  Cx , 

where  x  G  Mn,  y  G  Mp,  and  u  G  Mm  represent  the  state,  output,  and  control  input  respectively.  The  reference 
trajectory  to  be  tracked  r  G  W  is  generated  by  an  exosystem  of  the  form, 


(13) 


w  =  Sw 
r  =  Qw 

where  w  G  represents  the  exosystem  state  and  a  tracking  error  e  G  is  defined  as 

e  =  y  —  r  =  Cx  —  Qw 

In  what  follows,  we  consider  the  case  in  which  the  plant  model  (13)  is  over- actuated,  that  is  m  >  p. 


(14) 


(15) 
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Definition  III.l  (Tracking  Problem)  The  tracking  problem  for  (13)  and  (If.)  consists  in  finding  a  control 
law  such  that 

•  The  origin  of  (13)  is  asymptotically  stable  when  w(t)  =  0. 

•  For  any  initial  conditions  of  the  system  and  exosystem  the  tracking  error  (15)  converges  asymptotically 
to  0. 


The  full  information  tracking  problem  assumes  that  x(t)  and  w(t)  are  available  for  feedback  and  the 
control  input  u  is  given  by, 

u  =  —Kx  +  Lwy  (16) 

where  K  G  Mmxn  and  L  G  Wnxq.  The  full  information  tracking  problem  is  formally  defined  as  follows. 

Definition  III. 2  (Full  Information  Tracking  Problem)  Given  {A,  B,C,  S,Q}  find  K  G  Mmxn  and 
L  G  Wnxq  such  that 

•  (A-BK)  is  Hurwitz. 

•  For  any  initial  condition ,  the  trajectory  of 

x  =  ( A  —  BK)x  +  BLw  , 

w  =  Sw 

satisfies  lim^oo  || Cx(t)  —  Qw(t)\  \  =  0 

Theorem  III.l  below  due  to  Francis,  gives  the  conditions  for  the  solvability  of  the  problem  in  terms  of  the 
existence  of  a  solution  to  a  certain  linear  matrix  equation.15, 17 

Theorem  III.l  The  full  information  tracking  problem  is  solvable  iff: 

1.  (A,B)  is  stabilizable 

2.  3  matrices  II  G  M.nxq  and  T  G  Wnxq ,  solutions  of  the  Francis  equation 


nS  =  All  +  BY 
0  =  CU-Q. 


The  condition  that  (A,  B)  is  stabilizable  can  be  tested  using  standard  techniques.18  The  solvability  of 
the  Francis  equations  is  more  involved.  For  our  study,  we  will  make  use  of  Theorem  III. 2  which  gives  a 
sufficient  condition  for  its  solvability.16, 17 

Theorem  III. 2  If 

rank  ^  A  ^  j  =  n  +  p,  (19) 

for  any  X  G  spec(S),  then  the  Francis  equations  (18)  are  solvable  for  any  Q. 

The  solution  of  the  Francis  equation  defines  the  steady  state  of  the  system  on  which  the  tracking  error 
is  identically  zero,  and  the  input  required  to  maintain  the  system  in  steady  state.  This  steady  state  solution 
can  be  expressed  in  terms  of  solutions  of  the  Francis  equations  as 

xss  =  Uw,  uss  =  Tw.  (20) 

The  next  issue  we  need  to  address  is  to  design  a  controller  that  steers  the  trajectories  of  the  system  to  the 
steady  state.  Consider  the  following  change  of  coordinates 

x  =  x  —  xss,  u  =  u  —  uss.  (21) 
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By  substitution  into  equations  (13)  and  (15)  we  can  express  the  dynamics  and  tracking  error  in  the  new 
coordinates  as 

x  =  Ax  +  Bu 
e  =  Cx. 

Observe  that  by  construction  if  x(t)  =  0  then  the  state  and  input  is  equal  to  the  steady  state  solution  (20) 
and  the  tracking  error  is  identically  zero.  Thus,  to  solve  the  tracking  problem  we  need  to  stabilize  the  origin 
of  equation  (22).  This  is  achieved  by  means  of  state  feedback  choosing  a  feedback  gain  K  E  Mmxn  such  that 
(A  —  BK )  is  Hurwitz  and  by  applying  the  control  law 

u  =  —Kx.  (23) 

The  application  of  the  control  law  (23)  renders  the  steady  state  solution  xss  given  by  the  first  identity  in 
(20)  attractive,  while  the  conditions  of  Definition  III. 2  are  satisfied  by  construction.  Note  that  the  choice  of 
a  steady  state  control  is  completely  independent  from  the  choice  of  a  stabilizing  control  law.  This  separation 
is  key  for  the  methods  we  will  develop  for  steady  state  optimization  in  later  sections. 

Combining  the  steady  state  with  the  stabilizing  control,  we  can  express  the  control  law  as 


u  =  uss  +  u  =  Tw  —  K(x  —  lire)  =  —  Kx  +  (T  +  KlVjw.  (24) 

By  comparing  this  resulting  control  law  to  Definition  III. 2  we  observe  that  the  feedforward  gain  can  be 
written  as  L  =  (T  +  KT1)  . 


A.  Solving  the  Over- Actuated  Tracking  Problem 

It  is  shown  in  Ref.  15  that  if  Theorem  III. 2  holds  for  a  square  system  (m  =  p),  then  the  solution  of  the  Francis 
equations  is  unique.  If  the  system  is  over  actuated  (m  >  p )  then  the  solution  is  in  general  not  unique.  In  this 
case,  we  seek  to  find  all  solutions  to  the  Francis  equations  and  parameterize  them.  This  parameterization 
can  then  be  utilized  for  steady  state  optimization  and  satisfying  constraints  on  the  steady  state. 

We  start  rewriting  equations  (18)  as 


I  0 
0  0 


a2 


S 


x 


or 

R  =  AxX  -  A2XS.  (25) 

Equation  (25)  is  a  special  case  of  the  so-called  Hautus  matrix  equations.17  For  brevity  let  us  denote  the 
dimensions  of  the  matrices  X  and  R  as  fh  and  p,  respectively,  i.e.,  let 

fh  :=  dim(X)  =  (n  +  m)q,  p  :=  dim (R)  =  (n+p)q,  fh  >  p.  (26) 

We  introduce  the  mapping 

H  :  R™  i— >  RP,  H{X)  =  AxX  -  A2XS ,  (27) 

so  equation  (25)  can  represented  as  H(X)  =  R.  If  the  conditions  of  Theorem  III. 2  hold,  then  Ti  is  onto, 
which  implies  that  the  equation  H(X)  =  R  is  always  solvable  for  any  R  E  MP.  In  order  to  find  an  algebraic 
solution  to  the  equation  H(X)  =  R,  we  express  it  as  the  linear  system  of  equations 


Hx  =  r  (28) 

where  the  H  E  x  and  r  are  defined  as 

H  =  Iqxq  0  Ai  -  ST  0  A2,  x  =  col(X),  r  =  col (R),  (29) 

where  the  symbol  ®  denotes  the  Kronecker  product,  and  x  E  Mm  and  r  E  W  are  tall  vectors  constructed 
from  the  column  vectors  of  X  and  R  respectively.  The  matrix  H  El^xm  has  rank  p  by  assumption.  Every 
solution  of  equation  (28)  can  be  expressed  as 

x  =  xP  T  x°  (30) 
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where  xp  E  Mm  is  a  particular  solution,  and  x°  E  ker(iL).  Let  q  denote  the  dimension  of  the  kernel  of  iL,  i.e. 

q  :=  dim(ker (H))  =  dim(X)  —  dim (R)  =  fh  —  p.  (31) 

From  equation  (30)  it  is  evident  that  there  exists  an  infinite  number  of  solutions  to  equation  (28)  and  the 
dimension  of  the  solution  space  is  precisely  q.  This  observation  leads  us  to  Proposition  III.3. 

Proposition  III. 3  Assuming  that  H  E  Wxrn  in  (28)  has  rank  p,  then  there  exist  q  linearly  independent 
vectors  such  that  x°  in  (30)  can  be  expressed  as 

x°  =  ki#i  +  k2#2  +  •  •  •  +  k q6q,  where  9{  E  M,  i  =  1., . . . ,  q  (32) 

The  proof  of  Proposition  III.3  follows  elementary  arguments.  In  particular,  the  set  of  vectors  {k^}  is 
any  set  of  linearly  independent  vectors  spanning  ker(iL).  Using  equation  (32)  we  can  write  the  solution  of 
equation  (25)  as 

=  X  =  Xp  +  X°  (33) 

where  the  column  vectors  of  Xp  and  X°  can  be  found  from  xp  and  x°,  respectively.  Therefore,  analogously 
to  Proposition  III. 3,  there  exists  q  linearly  independent  matrices  Ab,  N2,  . ..,  Nq  E  M(n+m)x^  spanning 
ker(7Y)  such  that 

X°  =  N101  +  N202  +  •  •  •  +  NqOq.  (34) 

To  recap,  we  now  have  our  expression  for  all  matrices  II  and  T  that  satisfy  the  Francis  equations, 
parameterized  by  0  E  as  in  (34),  where  0  is  an  arbitrary  constant  vector.  Therefore,  with  a  slight  abuse 
of  notation  justified  by  the  linearity  of  the  parameter  #,  we  can  write  the  steady  state  solutions  as 

xss  =  Uw  =  IPw  + 11(6)  w  =  IPw  +  II  (w)0 

uss  =  Tw  =  Tvw  +  T^w  =  rp^  +  Y(w)0. 

The  control  law  can  be  expressed  as 

u  =  -Kx  +  (TP  +  T(<9)  +  K(W  +  H(0)))w. 

where  K  E  Mmxn  is  the  stabilizing  feedback,  and  0  E  is  an  arbitrary  constant  vector  that  defines  a 
particular  steady  state  solution. 


(35) 

(36) 


IV.  Optimization 

In  this  section,  we  employ  optimization  methods  to  choose  the  value  of  the  parameter  vector  9.  The 
optimization  objective  is  to  minimize  a  given  cost  function  on  a  moving  horizon  given  by  the  time  interval 
between  the  knot  points  of  a  reference  spline  trajectory,  generated  by  a  suitable  exosystem. 

The  cost  function  we  will  consider  depends  on  the  integral  of  the  steady  state  input  over  a  given  time 
interval  and  the  squared  2-norm  of  the  difference  between  9  and  the  optimal  choice  of  9  for  the  previous 
time  interval.  We  will  consider  both  an  offline  and  an  online  approach  to  solve  the  optimization  problem. 
When  using  an  offline  approach,  the  optimal  constant  9  is  computed  for  a  given  time  interval  during  the 
previous  one.  To  apply  the  optimal  #,  the  value  of  the  parameter  vector  9  is  simply  switched  at  the  start  of 
each  interval.  When  using  an  online  approach,  dynamics  are  assigned  to  9  such  that  during  each  given  time 
interval,  9  converges  continuously  to  the  offline  optimal  9. 

A.  Problem  Formulation 

The  system  that  we  will  consider  is  defined  in  Section  III  by  equations  (13)  through  (15),  with  the  assump¬ 
tion  of  over-actuation,  m  >  p,  and  that  the  full  information  tracking  problem  is  solvable.  The  following 
assumptions  define  the  class  of  reference  trajectories  generated  by  the  exosystem  (14). 

Assumption  IV.  1  The  exosystem  in  (14)  generates  cubic  splines  as  reference  trajectories  for  each  output 
of  the  system  in  (13). 
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Figure  2.  Controller  Structure  (the  ’’trigger”  is  a  trigger  signal  sent  at  t  =  tk  for  each  k  E  N) 


Assumption  IV. 2  There  exist  a  time  T  >  0  such  that  for  all  k  E  N  the  time  interval  [4,4+ 1)  between 
spline  knot  points  is  such  that  (tj. c+i  —  tk)  >  T. 


Assumption  IV.  1  is  fulfilled  using  a  chain  of  integrators  and  resetting  the  initial  conditions  at  the  knot  points 
of  the  splines.  The  purpose  of  Assumption  IV. 2  is  to  ensure  that  there  is  enough  time  to  complete  the  steady 
state  optimization  computations  for  the  next  time  interval.  By  choosing  T  large  enough,  we  also  ensure  that 
if  6  is  changed  at  the  beginning  of  a  time  interval,  the  stabilizer  will  have  sufficient  time  to  let  the  trajectories 
of  the  system  converge  to  the  steady  state  before  the  end  of  the  interval.  In  the  actual  implementation,  this 
will  be  achieved  by  tuning  T,  the  stabilizing  gain,  the  rate  of  convergence  of  6  to  the  offline  optimal  6  (for 
an  online  approach),  and  the  weights  in  the  cost  function. 

Let  6k  be  the  constant  steady  state  parameter  0  over  the  kth  time  interval  between  the  kth  and  (1  +  k)th 
reference  spline  knot  points,  i.e. 

6(t)  =  0k,  Vt  E  [4,4+ 1).  (37) 

We  will  consider  a  scalar  cost  function  Jk(6k )  :  i— ►  M  defined  over  the  kth  interval  between  the  reference 
spline  knot  points.  For  each  time  interval,  we  desire  to  find  the  value  such  that  the  given  cost  function 
Jk{6k )  is  minimized,  i.e. 

61  argmin  Jk{0k)  (38) 

Ok 

The  cost  function  we  will  consider  is  intended  to  penalize  the  input  energy  at  steady  state  over  the  given 
time  interval,  while  also  weighting  the  change  of  6k  in  squared  2-norm. 

The  cost  function  is  chosen  as 


J k  — 


«i  ftk+1 

-  / 

Atk  Jtk 


s(wk,ek,T)rruss(wk,ek,T)dT  +  -\6k  -  6+1 12 

A4 


(39) 


where  aq  E  M  and  ol^  E  M  are  positive  weighting  parameters,  Atk  =  (4+1  —  4)?  and  the  steady  state  input 
over  t  E  [4,4+i)  resulting  from  the  choice  6  =  6k  is  written  as 


uss(wk,9k,t)  =  (rp  +  T{6k))es(-t-t^wk, 


(40) 


where  wk  =  w(tk).  Obviously  by  increasing  the  ratio  o+ai  we  will  reduce  the  change  in  norm  between 
91  and  <+x.  However,  we  want  to  avoid  choosing  012/011  large  since  our  main  objective  is  to  minimize  the 
steady  state  input  energy. 

We  will  solve  the  optimization  problem  in  equation  (38)  with  an  offline  approach  and  an  online  approach 
for  both  unconstrained  and  constrained  optimization  to  enforce  constraints  on  the  input  at  steady  state. 
The  offline  approach  refers  to  changing  6  in  discrete  steps  as  the  exosystem  is  reset,  and  by  online  we  refer 
to  assigning  dynamics  to  6{t)  such  that  6(t)  converges  continuously  to  the  offline  optimal 

The  structure  we  seek  to  implement  is  depicted  in  Figure  2.  Notice  the  distinct  separation  between  the 
steady  state  control  or  design  and  the  stabilization,  this  results  from  the  specific  formulation  of  the  problem. 
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Figure  3.  Simplified  conceptual  example  for  offline  optimization  of  0. 


B.  Offline  Approach 

The  offline  approach  to  unconstrained  optimization  of  the  given  cost  function  Jk  involves  determining  the 
optimal  01  over  each  time  interval  of  the  reference  spline  and  applying  it  by  changing  0  in  a  discrete-time 
fashion.  Figure  3  illustrates  this  idea  and  also  points  out  that  switching  0  induces  discrete  changes  in  the 
steady  state  control,  and  thus  switching  in  uss(t )  will  induce  transients  in  u(t). 

To  achieve  unconstrained  offline  optimization  of  (39)  we  note  that  Jk  given  in  (39)  is  convex  in  uss ,  and 
uss  depends  linearly  on  6.  The  minimizer  0^  is  found  solving  V#fc  Jfc(0£)  =  0,  where  WokJk  denotes  the 
gradient  of  Jk  with  respect  to  Ok-  To  see  why  this  is  the  case,  we  compute  the  gradient, 

VekJk  =  Lk  +  Rk@k,  (41) 


where 


Lk  =  fk  ftj1  ^MVpwdr  -  fkek~v  Rk  =  tii?+lr»rH^+t4 


Atk  Jtk 


(42) 


To  solve  Jk{9%)  =  0>  we  simply  need  to  find  the  inverse  of  Rk,  and  obtain 


01  = 


-R^Lk. 


(43) 


Since  ot2  >  0,  and  Rk  is  formed  by  ot2  times  the  identity  matrix,  plus  a  positive  semi-definite  matrix  resulting 
from  the  integration  of  a  matrix  times  its  transpose,  then  Rk  is  positive  definite.  Hence,  R^1  is  guaranteed 
to  exist.  Furthermore,  since  Rk  >  0  then  the  Hessian  of  Jk(0j. c)  with  respect  to  Ok  is  positive  definite,  thus 
0^  is  the  global  minimizer  of  the  cost  function.19  We  therefore  have  a  valid  solution  to  the  unconstrained 
offline  optimization  problem  and  the  solution  is  expressed  by  equation  (43). 


C.  Online  Approach 

Motivated  by  the  need  to  reduce  the  transient  behavior  that  can  occur  in  the  control  input  and  the  plant 
output  as  a  result  of  switching,  we  desire  to  take  an  approach  to  change  0  continuously.  We  solve  the 
problem  by  assigning  dynamics  to  0k(t)  such  that  0k(t)  rapidly  converges  to  the  offline  optimal  0%.  Consider 
the  gradient-based  update  law 

Ok  =  -1^7  ejk,  te[tk,tk+ 1),  (44) 

where  7  >  0  is  a  design  parameter.  Recall  that  the  cost  function  Jk  changes  for  each  time  interval,  which 
makes  the  gradient-based  update  law  change  discretely  at  the  time  of  each  reference  spline  knot  point.  We 
will  show  that  this  choice  of  0k(t)  dynamics  renders  0 £  an  asymptotically  stable  equilibrium,  while  the  gain 
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V'SS 


Figure  4.  Simplified  conceptual  example  for  online  optimization  of  0. 


7  affects  the  convergence  rate  of  6k(t)  to  0 To  show  this  we  start  by  defining  the  following  change  of 
coordinates 

Ok  =  Ok  -01  (45) 

Recall  that  by  definition  {Lk  +  RkO^)  =  0  and  therefore  the  dynamics  in  the  new  coordinates  can  be  expressed 
as 

Ok  =  —^RkOk-  (46) 

Note  that  since  Rk  >  0,  then  —'yRk  has  all  negative  eigenvalues  and  6  =  0  is  an  asymptotically  stable 
equilibrium.  Furthermore,  7  scales  the  eigenvalues  of  —7/4  and  thus  the  rate  of  convergence  of  Ok  to  0 1  is 
proportional  to  7. 

Figure  4  shows  in  a  simplified  way  how  using  a  continuous  0(t)  will  produce  a  continuous  uss(t ),  which 
will  help  reduce  transients  in  u(t). 

V.  Constrained  Optimization 

It  is  well  known  that  in  most  control  applications  there  are  constraints  on  the  states  and  inputs  of  the 
system  to  be  controlled  and  the  air-breathing  hypersonic  vehicle  is  certainly  no  exception.  In  this  section,  we 
will  add  constraints  to  the  steady  state  input  for  the  system  we  are  considering.  In  particular,  we  formulate 
our  criteria  as  follows: 

•  Satisfy  possibly  asymmetric  constraints  on  some  or  all  of  the  inputs  at  steady  state  (uf) 

•  Minimize  the  energy  of  the  input  at  steady  state,  as  a  secondary  objective. 

We  will  assume  that  it  is  possible  to  satisfy  this  criteria.  That  assumption  is  made  for  illustrative  purposes 
that  the  constrained  problem  fits  directly  into  our  framework.  To  simplify  the  problem,  we  will  consider  only 
a  constraint  on  one  of  the  inputs,  the  ith  input.  An  offline  and  an  online  approach  will  be  considered.  Since 
the  online  approach  utilizes  the  results  of  the  offline  approach,  the  online  approach  is  more  aptly  referred  to 
as  an  online  extension. 

First,  we  need  to  set  up  some  assumptions,  notation  and  terminology.  We  define  constraints  on  the  ith 
input  at  steady  state  as  follows:  given  scalars  u  <  u,  we  desire  a  control  at  steady  state  uss  such  that 

uss  G  U ,  where  U  =  {uss  G  Mm  :  u  <  uf  <  u},  (47) 

where  uf  denotes  the  ith  steady  state  input.  We  define  0k(t)  to  be  a  feasible  parameter  trajectory  at  time 
t  G  [tfc,tfc_|_  1)  if  u  <  uf(wk,0kR)  <  u.  Furthermore,  let  us  define  the  set  ©&  as 

Qk  =  {6  :  u  <  uf  <  u,  uss  =  (Tp  +  T(0))es{-t~tk)wk ,  Vi  e  [4,4+ 1]-}  (48) 
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and  refer  to  Ok  as  the  feasible  parameter  set  over  the  time  interval  [tk,tk+ 1).  Note  that  even  though  U  is 
not  time  varying,  the  boundary  dOk  of  the  feasible  set  Ok  does  in  general  vary  between  time  intervals.  Also 
note  that  by  continuity  Ok  G  Ok  implies  that  for  some  t£  >  0,  0k(t)  is  feasible  over  \tk+i,tk+i  + 1£).  We  will 
assume  that  Assumptions  IV.  1  and  IV. 2  still  hold.  In  addition,  we  assume  that  there  is  sufficient  control 
authority  to  solve  the  constrained  optimization  problem: 

Assumption  V.l  Given  the  full  information  tracking  problem  as  described  in  Section  III )  Vfc  G  N,  3  Ok  G  Ok 
and  a  scalar  e  >  0  such  that  (u  +  e)  <  uf(0k,t)  <  (u  —  e),  V  t  G  [£&,  tk+ 1]. 

The  constrained  problem  can  now  be  expressed  as  finding  the  optimal  0 1  over  the  kth  time  interval  such 
that 

01  =  arg  min  Jk(0k)  (49) 

where  Ok  is  defined  by  (48)  and  Jk(0k)  is  defined  by  (39).  As  with  the  unconstrained  case,  we  will  consider 
both  an  online  and  an  offline  approach  to  solving  this  problem.  We  will  change  0  in  discrete  steps  for  the 
offline  case  and  assign  the  dynamics  of  0k  for  the  online  case  such  that  0k(t)  converges  to  0 1  without  becoming 
unfeasible  at  any  point  in  time. 

A.  Offline  Approach 

The  problem  of  finding  an  offline  constrained  optimal  0k  over  the  time  interval  \bk,tk+ 1]  will  be  split  into 
two  steps: 

•  Step  1:  Find  an  initial  guess  0®  G  Ok  for 

•  Step  2:  With  such  an  initial  guess  for  the  optimization,  find  a  0 ^  that  minimizes  the  cost  function  Jk 
defined  by  (39)  in  the  sense  of  (49),  i.e.  optimize  such  that  01  G  Ok- 

Since  our  main  objective  is  to  ensure  that  a  given  Ok  is  in  the  feasible  set  ©&,  we  must  have  a  way  to 
test  for  such  feasibility.  Let  us  define  a  2nd  order  polynomial  of  uf  that  is  negative  iff  uf  is  within  the 
constraints,  i.e. 

c(uf )  =  uf2  +  C\uf  +  Co,  such  that  c(u)  =  c(u)  =  0,  and  c(uf )  <  0  \/uf  G  {uf  :  u<  uf  <  u}. 

(50) 

Since  w(t)  is  a  polynomial  in  t  then  uf  is  a  polynomial  in  t  which  allows  us  to  also  express  c(uf )  as  a 
polynomial  in  t  i.e.  c{rk)  with  a  slight  abuse  of  notation,  where  Tk  =  t  —  tk-  If  the  roots  of  c(rk)  are  within 
[0,  At^],  then  the  constraints  will  be  violated  within  the  kth  time  interval.  We  have  thus  come  up  with  a 
feasibility  check  that  can  also  tell  us  at  what  point,  or  points,  in  time  an  unfeasible  Ok  allows  the  steady 
state  input  to  violate  a  constraint. 

Our  next  task  is  to  find  a  0  in  Ok-  We  will  accomplish  this  by  using  a  gradient  optimization  technique. 
The  method  involves  setting  up  a  cost  function  that  makes  the  cost  small  if  the  input  is  within  the  constraints 
but  rapidly  gets  larger  if  the  input  is  near  or  outside  the  constraints.  Since  by  Assumption  V.l  there  exists 
an  input  that  is  feasible  over  the  given  time  interval,  then  there  should  exist  an  input  that  makes  this  cost 
function  small.  First  a  cost  function  is  defined  as  follows 

i  rtk+i 

J(t>  =  ^vT  a3<p(uf)dT  (51) 

Atk  Jtk 

where  0^3  >  0  is  a  weighting  gain.  We  choose  </>(•)  as  a  convex,  positive  definite  function  with  <j>(u)  =  <j)(u). 
We  will  use 

O/l/SS  _  n  i\ 

4>(uf)  =  l{uf)2u ,  0  <  V  G  N,  where  l{uf)  =  •  (52) 

Basically  we  construct  a  line,  /(•),  that  intersects  -1  and  1  at  u  and  u  respectively,  then  </>(•)  is  /(•)  raised 
to  an  even  power  to  create  the  desired  convex  function  that  is  exactly  1  at  u  and  u.  Given  Assumption  V.l, 
choosing  a  large  enough  v  guarantees  that  minimizing  J ^  with  respect  to  0k  will  result  in  a  ^  in  6/;.  Let  t^ 
be  defined  as 

4  =  min({te,tfc+i})  (53) 

where  te  is  such  Ok  is  feasible  over  [tk,te)  and  unfeasible  at  te  if  Ok  is  unfeasible  at  any  time  >  £&,  else  te  =  oo. 
Then  to  get  a  Ok  that  is  in  Ok  we  only  need  to  iteratively  minimize  J f  until  t^  =  tk+\.  This  is  method  is 
described  by  Algorithm  V.2. 
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Algorithm  V.2  (Step  1)  Let  0(_t  G  @/c-i  denote  the  minimizer  for  the  previous  reference  spline  interval. 
By  continuity  0(_%  is  feasible  over  at  least  [tkHk  + t£)  for  some  t£  >  0. 

1.  Assign  6®  =  0(_1.  If  0^  G  ©&,  do  nothing  and  return  6®  =  Ok- 1- 

2.  Otherwise,  iteratively  minimize  J $  until  we  get  t^  =  tk+i,  he.  we  have  a  6®  that  is  in  ©&.  Return  6®. 

The  next  thing  we  need  to  do  is  optimize  0  over  [4,h+i)  without  violating  the  constraints  on  uf*.  We 
will  do  this  with  a  gradient  method  initialized  with  0 £  employing  the  same  cost  function  as  we  used  for 
the  unconstrained  steady  state  optimization,  i.e.  Jk  from  equation  (39).  To  satisfy  the  constraints  in  an 
optimization  algorithm  we  implement  Algorithm  V.3. 

Algorithm  V.3  (Step  2)  Let  6®  be  the  initial  guess  for  0(,  obtained  using  Algorithm  V.2. 

Iteratively  optimize  Jk  (39)  (gradient  method)  until  we 

1.  Have  minimized  Jk  to  within  tolerance  limits,  stop  iterating  and  return  0(. 

2.  Taking  a  step  in  the  gradient  search  that  is  smaller  in  norm  then  a  given  threshold  causes  0 (  not  to  be 
in  &k,  stop  iterating  and  return  the  feasible  0(  without  further  gradient  search  steps. 

B.  Online  Extension 

When  dealing  with  the  online  constrained  optimization,  we  take  a  different  approach  then  the  one  taken  for 
the  unconstrained  case.  We  call  this  an  extension  because  we  will  still  use  the  offline  optimal  algorithm  to 
find  a  reference  for  Ok-  For  the  constrained  online  optimization,  we  will  follow  the  following  steps: 

•  We  will  assume  that  we  have  used  the  offline  algorithm  and  found  0 (  during  the  previous  interval 
\fk  —  1  1  tk) 

•  Assigning  0k(t)  to  steer  (Ok  —  0(f)  to  zero,  arbitrarily  fast,  without  letting  0k(t)  become  unfeasible  at 
any  time. 

Recall  that  0(_1  is  feasible  at  the  next  interval  for  at  least  some  [tkHk  +  te],  t£  >  0.  Therefore,  we  have 
at  least  te  units  of  time  to  bring  0k(t)  to  a  small  neighborhood  of  0(  before  it  becomes  unfeasible  (tifs(£) 
violates  its  constraints). 

Similarly  to  the  online  case  in  the  unconstrained  problem,  we  will  be  considering  assigning  a  dynamics 
to  0k(t)  such  that  0k(t)  converges  to  0(  faster  then  the  dynamics  of  the  system  .  This  assumption  is  made 
in  order  to  justify  that  the  changes  in  0  will  only  cause  transients  in  the  state  and  input  that  the  stabilizing 
controller  is  able  to  render  negligible  within  [tk,tk+ 1). 

To  remedy  the  problem  of  possible  infeasibility  during  transients  in  0k(t)  when  using  a  simple  feedback 
law  in  Ok  (t)  as  was  done  for  the  unconstrained  online  optimization,  we  suggest  adding  a  barrier-type  function 
to  the  dynamics  of  Ok .  Details  are  skipped  for  reasons  of  space,  and  will  be  given  in  the  final  version  of  the 
paper. 


VI.  Simulation  Results 

The  control  and  optimization  methods  described  in  the  previous  chapters  were  applied  to  the  hypersonic 
air-breathing  vehicle  model  described  in  the  Section  I.  The  results  presented  in  this  section  are  mostly 
preliminary.  Further  tuning  might  potentially  improve  the  results  significantly.  The  physical  meaning  of  the 
inputs  has  been  briefly  discussed  in  the  introduction  of  the  model  but  here  we  need  to  consider  them  in  slightly 
more  detail  to  make  our  optimization  and  assigning  of  constraints  meaningful.  The  diffuser  area  ratio  A&,  can 
obviously  only  vary  between  [0, 1].  We  chose  to  use  a  more  conservative  range,  [0.14,  0.735].  The  range  for 
the  total  temperature  change  across  the  combustor  will  be  considered  to  be  [500  degR,  3900  degR].  Finally, 
the  elevator  deflection  should  remain  within  the  range  [—30°,  30°].  The  addition  of  actuator  dynamics  to 
the  model  takes  care  of  the  bandwidth  limitation  of  the  inputs.  Because  of  the  large  differences  in  input 
ranges,  we  scale  the  B  matrix  to  normalize  their  relative  effect  of  the  input  on  the  optimization  methods. 
The  scaling  factors  and  parameters  used  for  the  steady  state  optimization  methods  are  presented  in  Table 
3,  and  Table  2  shows  the  trim  state  used  for  the  linearization  of  the  model  dynamics. 

We  considered  four  cases  in  our  simulations: 
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Figure  5.  Hypersonic  air-breathing  vehicle  model:  Reference  for  velocity  and  altitude. 


•  Offline  unconstrained  steady  state  optimization  applied  on  the  linear  plant  (OFF-NC-L) 

•  Offline  constrained  steady  state  optimization  applied  on  the  linear  plant  (OFF-C-L) 

•  Online  constrained  steady  state  optimization  applied  on  the  linear  plant  (ON-C-L) 

•  Online  constrained  steady  state  optimization  on  the  nonlinear  plant  (ON-C-NL) 

We  applied  the  offline  unconstrained  optimization  to  help  tune  the  stabilizer  and  get  fair  results  while  allowing 
Ad  to  violate  its  constraints.  Then  we  considered  the  offline  constrained  optimization  to  prevent  violating 
the  constraints  on  Ad.  The  online  constrained  optimization  was  used  to  reduce  transients  in  the  response 
observed  for  the  offline  case.  Finally,  the  online  constrained  optimization  was  applied  to  the  nonlinear  model. 

For  the  stabilizing  controller,  we  chose  to  use  state  feedback  with  a  gain  K  as  in  the  examples  given  in 
the  previous  chapters.  To  tune  the  feedback  gain  we  use  LQR  methods  on  the  linear  system  with  the  scaled 
B  matrix  and  the  following  weighting  matrices, 

R  =  hx3 

Q  =  diag(10.0,  1(T3,  1(T3 *,  1.0,  1(T3,  1(T6,  1(T6,  1(T6,  1(T6,  0.0,  0.0,  0.0). 

For  the  unconstrained  cases: 
a  i  =  l  0^2=0. 1  7  =  100 

For  the  constrained  cases: 

o,i  =  1  <a2  =  10  7  =  0.3  u  =  -0.21  u  =  0.385  constrain  uf 

a3  =  1  v  =  4 
Input  scaling  factors: 

5e:  1/10  A  T0:  1/1000  Ad:  1/0.7 

Table  3.  Parameters  for  the  optimization  and  scaling  of  the  inputs.  The  constraints  are  presented  w.r.t. 

unsealed  input  as  deviations  from  the  trim  condition  and  are  conservative. 


The  reference  velocity  and  altitude  that  were  chosen  are  shown  on  Figure  5.  We  chose  an  approximation 
of  a  step  of  1000  ft/sec  in  velocity  and  10000  ft  in  altitude  over  time  span  of  93  s  (within  1%  of  the  initial 

and  final  value).  Those  are  compatible  references  for  velocity  and  altitude  to  the  ones  used  in  Ref. 6.  Figure 

6  shows  the  tracking  error,  which  is  small  with  respect  to  the  amplitude  of  the  reference.  The  velocity 
and  altitude  references  are  of  the  order  of  103  and  104,  respectively,  while  the  error  is  only  10“ 1 *  and  10° 
respectively.  As  expected,  the  error  for  the  nonlinear  model  is  greater  then  for  the  linear,  but  it  is  still  very 
small. 

Figure  7  shows  uss.  We  observe  that  for  the  unconstrained  case  we  grossly  violate  the  constraints  on  Ad . 
However,  we  also  see  that  all  the  constrained  cases  do  not  allow  Ad  to  become  unfeasible.  The  variations  in 
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Figure  6.  Tracking  error  for  velocity  and  altitude. 


Se  are  also  well  within  reasonable  amplitudes.  ATq  becomes  quite  large  but  that  can  possibly  be  adjusted 
by  further  tuning.  Figure  8  shows  the  inputs  to  the  model  of  the  hypersonic  air-breathing  vehicle,  notice  the 
transients  caused  by  changing  9.  As  expected,  the  online  approach  greatly  reduces  the  transients.  In  fact 
the  transients  are  reduced  enough  that  for  the  online  cases  the  transients  in  A d  not  to  exceed  its  constraints. 

Figure  9  shows  9  for  the  cases  we  considered.  We  have  an  eight  dimensional  9  parameter  so  visualizing 
its  effects  directly  from  the  graph  is  difficult  at  best.  Because  of  this  difficulty,  the  only  comments  we  will 
make  on  this  graph  are  that  we  observe  that  the  different  methods  we  use  to  choose  9  give  different  results, 
and  that  the  online  methods  give  a  continuous  #,  while  the  offline  methods  change  9  discretely.  The  rapid 
convergence  of  6k(t)  for  the  online  cases  to  the  offline  optimizer  91  makes  it  difficult  to  distinguish  between 
the  two  approaches  on  Figure  9. 

Finally,  we  plot  some  of  the  states  that  are  not  outputs  on  Figure  10.  The  angle  of  attack  is  of  great 
importance  since  the  vehicle  may  not  function  if  the  amplitude  of  a  becomes  too  large,  i.e.  goes  outside 
[—5°,  5°].  Figure  10  shows  that  a  remains  within  about  1.5°  from  the  trim  condition  for  the  constrained 
cases.  It  does  seem  to  make  rather  abrupt  changes  but  note  that  the  time  scale  is  considerably  long  with 
respect  to  the  dynamics  of  the  system.  The  pitch  angle  it  varies  only  up  to  about  2-3°  from  the  trim  which 
also  is  not  great  enough  to  cause  concern.  Observe  also  that  the  flight  path  angle  (FPA)  is  reasonable 
and  evolves  without  making  abrupt  changes  in  amplitude.  Finally,  the  generalized  flexible  modes  vary  in 
amplitude  no  more  then  what  was  observed  in  Ref. 6. 

VII.  Conclusions 

The  model  of  the  air-breathing  hypersonic  vehicle  presents  a  significant  control  challenge  as  was  discussed 
in  this  study.  Linearization  was  used  to  produce  a  simplification  the  model  which  yields  itself  better  to  control 
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(b)  AT0 


(c)  Ad 


Figure  8.  Inputs  to  the  model  of  the  hypersonic  air-breathing  vehicle  ( uv ). 
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Figure  9.  The  parameter  6  for  the  cases  under  consideration. 


design.  A  control  method  was  presented  that  allows  designing  a  steady  state  and  a  stabilizing  controller 
separately  for  tracking  control  of  this  overactuated  system.  Using  state  feedback  for  stabilizing  control  and 
optimization  methods  for  steady  state  control,  both  unconstrained  and  constrained,  we  were  able  to  produce 
promising  results  in  simulation.  Those  results  show  very  tight  tracking  performance  for  reference  trajectories 
that  are  more  demanding  that  those  presented  in  Ref.  20.  Also  we  showed  that  we  were  able  to  produce  near 
perfect  tracking  under  steady  state  input  constraints  even  when  applying  the  controllers  to  the  full  nonlinear 
model  of  the  air-breathing  hypersonic  vehicle. 

The  work  presented  in  this  study  opens  up  a  door  to  further  research  on  methods  using  the  regulator 
approach  to  the  tracking  problem  for  over-actuated  systems.  Robustness  analysis  of  the  control  methods, 
and  a  detailed  analysis  of  the  stability  and  rate  of  convergence  for  the  constrained  optimization  methods 
are  examples  of  potential  research  topics.  We  should  also  note  that  in  this  paper,  we  do  not  necessarily 
assume  we  have  a  differentially  flat  system,  which  is  a  standard  assumption  made  for  conventional  trajectory 
optimization  methods  applied  to  over-actuated  systems.11,13  Finally,  detailed  geometric  analysis  of  the 
regulator  problem  formulation  and  solution  may  lead  to  analogous  approaches  to  non-linear  control,  then  a 
new  non-linear  control  method  could  potentially  be  developed. 
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